arXiv:1503.07505vl [physics.bio-ph] 25 Mar 2015 


AN IMMERSED BOUNDARY MODEL OF THE COCHLEA WITH 
PARAMETRIC FORCING 

WILLIAM KOt AND JOHN M. STOCKIEt 


Abstract. The cochlea or inner ear has a remarkable ability to amplify sound signals. This is understood 
to derive at least in part from some active process that magnifies vibrations of the basilar membrane (BM) and 
the cochlear partition in which it is embedded, to the extent that it overcomes the effect of viscous damping from 
the surrounding cochlear fluid. Many authors have associated this amplification ability to some type of mechanical 
resonance within the cochlea, however there is still no consensus regarding the precise cause of amplification. Our 
work is inspired by experiments showing that the outer hair cells within the cochlear partition change their lengths 
when stimulated, which can in turn cause periodic distortions of the BM and other structures in the cochlea. 
This paper investigates a novel fluid-mechanical resonance mechanism that derives from hydrodynamic interactions 
between an oscillating BM and the surrounding cochlear fluid. We present a model of the cochlea based on the 
immersed boundary method, in which a small-amplitude periodic internal forcing due to outer hair cells can induce 
parametric resonance. A Floquet stability analysis of the linearized equations demonstrates the existence of resonant 
(unstable) solutions within the range of physical parameters corresponding to the human auditory system. Numerical 
simulations of the immersed boundary equations support the analytical results and clearly demonstrate the existence 
of resonant solution modes. These results are then used to illustrate the influence of parametric resonance on wave 
propagation along the BM and explicit comparisons are drawn with results from another two-dimensional cochlea 
model. 


Key words, cochlea, basilar membrane, immersed boundary method, fluid-structure interaction, parametric 
resonance. 
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1. Introduction. The cochlea or inner ear is a fluid-filled, spiral-wound cavity that is the 
central source of frequency selectivity in the hearing system of humans and many other animals. 
The cochlea is divided on its primary axis into two main fluid chambers, the scala tympani and 
scala vestibuli, by a structure known as the cochlear partition (CP). The CP itself consists of an 
elastic membrane known as the basilar membrane (BM), on top of which is mounted the organ of 
Corti containing the mechanosensitive hair cells that are the primary sensory receptors in the ear. 
Figure [TTlI contains a picture of an unwound cochlea with a simplified view of the CP showing only 
the BM. Sound vibrations entering the outer ear are transferred to the cochlear duct and BM by 
the stapes, and then propagate from base to apex along the BM and the surrounding fluid. The 
BM has a fine frequency tuning ability that distinguishes between different frequencies by localizing 
the peak amplitude of incident travelling waves. This “place theory” was developed and validated 
experimentally by von Bekesy [M] and shows that the peak location is closer to the base for high 
frequencies and to the apex for low frequencies, as illustrated in Figure [ITT] 

A second defining feature of the cochlea, which sets it apart from other sensory receptors, is its abil¬ 
ity to amplify extremely weak sound signals that would otherwise be immediately damped out due 
to viscous dissipation in the cochlear fluid. Many authors have attributed this amplification ability 
to some active process related to resonance, which experiments have connected with mechanical 
properties of various structures making up the CP In particular, the outer hair 

cells in the organ of Corti are stimulated by BM deflections caused by pressure waves travelling 
through the cochlear fluid. The hair cell stimulation leads to either somatic motility , wherein the 
hair cell changes its length in response to electrical signals induced when the hair bundle on its tip 
is deflected; or active hair bundle motility in which the hair bundle itself generates additional forces 
that initiate a shearing action between the BM and the overlying tectorial membrane. Both of 
these effects are believed to contribute to the cochlear active process, but there is still no consensus 
in the literature regarding the precise cause of amplification 
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Fig. 1.1. Diagram of an uncoiled cochlea and BM, showing the relative location of maximal BM response for 
several sound frequencies in the human audible range. This is a simplified view of the cochlear partition, depicting 
only the BM and omitting other structural components such as the organ of Corti, outer hair cells, and tectorial 
membrane. Source: \25i Fig. 2] (distributed under the Creative Commons License). 


Many mathematical and computational models of the cochlea have appeared since the seminal 
work of von Bekesy [53]. The earliest two-dimensional models of the cochlea describe the BM 
as a collection of damped mass-spring systems and also reduce the fluid dynamics to a simple 
linear potential flow [u na uni- The spring constant decreases exponentially along the BM from 
base to apex, which coincides with BM stiffness values measured in experiments [54]. The BM is 
treated as a passive structure to which is applied a sinusoidal external forcing term that mimics the 
input of sound energy from the stapes. These models give predictions of BM dynamics that agree 
qualitatively with the behaviours observed by von Bekesy. Inselberg and Chadwick [2T] proposed 
a similar model in which the BM is represented as an Euler-Bernoulli beam, and showed not only 
that the place principle still holds but also that fluid viscosity is required to obtain travelling wave 
solutions along the BM as opposed to simply standing waves [7:. Pozrikidis |47j revisited this last 
approach by replacing the sinusoidal stapes motion with a point source at the stapes position and a 
point sink at the round window, and then solving the resulting equations using a boundary integral 
method. Another noteworthy class of models based on transmission line equations was introduced 
in the pioneering work of Zweig [SB] and de Boer [12] and has since been applied in many more 
recent studies such as mm- 

To obtain a more realistic model of the fluid dynamics in the cochlea as well as the hydrodynamic 
interactions between the fluid and BM, several authors have exploited the immersed boundary (or 
IB) method. The IB method was originally developed by Peskin to simulate blood flow in the 
beating heart [15] and has since been applied to many problems in biofluid dynamics, including the 
cochlea. LeVeque et al. .30, 31] employed an IB model in which the fluid obeys the linearized Stokes 
equations and the elasticity of the cochlea is treated using simple elastic springs. They derived an 
asymptotic solution for travelling waves along the BM, from which they drew conclusions regarding 
the effects of fluid viscosity on these waves. Beyer and LeVeque [1] performed numerical simulations 
of a related IB model, with the primary difference being that their fluid obeys the full (nonlinear) 
Navier-Stokes equations. All of the aforementioned immersed boundary models approximate the 
geometry of the cochlea and BM by a straight (uncoiled) configuration. Although the curvature of 
the cochlear duct has a relatively small influence on the fluid dynamics m, there is nonetheless 
some evidence to suggest that within the most tightly coiled apical region of the BM that is 
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stimulated by the lowest frequency sounds, curvature cannot be ignored [34] . To this end, a much 
more detailed IB model capturing the full 3D geometry of the cochlea was developed in m that 
reproduced important features of BM dynamics. 

The goal of this paper is to use the IB method to investigate possible resonant phenomena that 
contribute to mechanical amplification of basilar membrane oscillations in the cochlea. Our study 
is inspired by three main observations: first, direct experimental evidence that outer hair cells 
embedded within the organ of Corti in the CP undergo periodic contractions when the ear is 
stimulated by sound waves [501 [531 EE]; second, the suggestion by several authors that these hair 
cell contractions can in turn modulate the stiffness in the BM m E3 [41] which we take as 
an assumption; and third, the analytical results of Cortez et al. m that uncovered parametric 
instabilities arising from the fluid-structure interaction in internally-forced immersed boundaries 
with a time-varying stiffness parameter. Taken together, these three observations suggest that 
there is merit in investigating the hypothesis that internal forcing in the basilar membrane due to 
sound stimulation can engender parametric resonance in an IB model of the cochlea. 

One of the main contributions of this work is to extend the previous parametric resonance analysis 
for an elastic membrane with a periodically-varying stiffness parameter m to the case where 
stiffness depends on both time and space. Our results may also contribute new insights into 
the understanding of fine-frequency tuning and amplification in the cochlea since many previous 
cochlea models (described above) treat the outer hair cell contractions as an external forcing, 
over-simplify the fluid mechanics, and/or ignore the fluid-structure interaction. For example, it is 
common for papers to question the ability of linked mechanical oscillator models to capture cochlear 
amplification by arguing that damping due to viscosity of the cochlear fluid is simply too large EH- 
This argument is valid if the active process driving the cochlea appears as an external forcing in the 
model, but not when the force acts internally through a parameter - in the latter case, parametric 
resonance can occur in which large amplitude oscillations (unbounded in the idealized linear case) 
arise even in the presence of damping [8]. Finally, we draw a distinction between our approach and 
other models incorporating resonance effects that act locally, such as the transmission line models 
based on coupled mechanical oscillators where each line element resonates individually. In contrast, 
our model exhibits global resonance owing to the non-local nature of the coupling between the BM 
and the surrounding fluid. 

2. Mathematical Model. Following the IB cochlea model derived by Peskin in [44l and 
developed in more detail in |31 . we consider a simple 2D geometry pictured in Figure [5TT1 in which 
the cochlear duct of length L is treated as a rectangular strip fl = [0 ,L\ x M, along the center of 
which lies the BM. This work is primarily concerned with the effects of parametric forcing on BM 

(+oo) 



Fig. 2.1. Geometry of the 2D immersed boundary model for the cochlea. The deformed BM is represented by 
a solid blue line and the flat equilibrium state by a red dotted line. To simplify the analysis, even symmetry is 
imposed across x = 0 and x = L so that the solution is defined for all a;Gl, although only values of x £ [0, L] have 
physical relevance. 


oscillations and so we simplify the model by isolating the membrane from any boundary effects 
due to the cochlear walls. For the purposes of the mathematical analysis, we take the depth of the 
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cochlear duct to be infinite (y —»• ±oo) which is consistent with the short-wave approximation [511 
and has been used previously by various authors USUI!- This assumption neglects effects such 
as coherent backscattering localized near the travelling wave peak IE!, but we have found that 
this does not significantly impact on the results of our stability analysis. Furthermore, we impose 
Neumann boundary conditions in the x-direction which is an assumption that has also been justified 
for other cochlea models mm- 

The fluid in which the BM is immersed is governed by the incompressible Navier-Stokes equations 


( 2 . 1 ) 

( 2 . 2 ) 


<9u 

~dt 


u • Vu 


-Vp + /rAu + f, 


V • u = 0, 


where u(x, t ) is velocity, p(x, t) is pressure, p is density, and p, is viscosity. The elastic force exerted 
on the fluid by the membrane is given by 

(2.3) f (x,i)=/ K(s, t) (X 0 — X) (5(x — X) ds, 

Jo 

where <5(x) is the two-dimensional Dirac delta function, X(s,f) is the position of the membrane 
parameterized by the Lagrangian coordinate s G [0,L], and Xo(s) = (s, 0) is the horizontal equi¬ 
librium or rest state. This force can be thought of as arising from a membrane that is connected 
to its resting position via a series of Hookean springs with “spring constant” or stiffness K(s,t). 
The elastic stiffness parameter is a function of time and space given by 

(2.4) K ( s,t ) = ere _As (l + 2 r sin (cot )), 

where a is the time-averaged elastic stiffness constant (units of gcm -2 s -2 ) and A captures the 
spatial variation in stiffness along the BM. The value of A ss 1.4 cm -1 has been determined for 
a human cochlea experimentally by von Bekesy |54i and others, based on the observation that 
the BM stiffness at the apex is approximately two orders of magnitude smaller than that at the 
base, and that it decays roughly exponentially. The time-dependent factor in the stiffness encap¬ 
sulates the parametric forcing with amplitude r and frequency to arising from outer hair cells that 
contract/expand in response to BM oscillations [fU [HI ED]. Note in particular that the forcing 
frequency ui is taken to be the same as that of the input sound signal and is also constant in space, 
which assumes that the outer hair cells contract in synchrony along the entire CP. This is in con¬ 
trast with most other models of the cochlea that consider hair cell contractions in response to local 
stimuli, which would correspond to a stiffness parameter having spatiotemporal dependence that 
is not separable. There is nothing in our analytical approach that prevents us from considering 
this situation, but we nonetheless restrict our attention to the separable form in (12.41) for the sake 
of simplicity. 

We close the system of equations with an expression for the membrane motion 

<9X f 

(2.5) — =u(X,t) = J u(x,t)5(x —X) dx, 

which is equivalent to the fluid no-slip condition along the immersed boundary. 

To simplify the analysis in the remainder of this paper, we non-dimensionalize the problem by 
introducing the following scalings 

Lx t „ _ LX Us 

(2.6) x = —, t=— , u = U c u, p = P cP , X=-, s = —, 

7T U) 7T 7T 

where the tildes denote a non-dimensional quantity. The characteristic velocity and pressure scales, 
U c and P c , will be specified shortly. The horizontal extent of the rescaled domain D = [0,7r] x E 
is chosen for reasons of mathematical convenience, in order to eliminate a factor of 7r that would 
otherwise appear in the solutions derived in section [3] 
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Substituting the above variables into the governing equations m m yields 


(2.7) 

<9u 

+ u • Vu = — Vp + j/Au + f, 
dt 

( 2 . 8 ) 

V • u = 0, 

(2.9) 

f(*,?)=/ A(X 0 -X)d(X 
J 0 

( 2 . 10 ) 

A'(s,f) = Ke -as (l + 2rsint), 

( 2 . 11 ) 

A =fl( x,o. 

dt 


The velocity and pressure scales are chosen as 

( 2 . 12 ) = ^ and P c = ^, 

7T 7T Z 


to reduce the number of dimensionless parameters appearing in the equations to four, consisting 
of r plus the three new quantities 


(2.13) 


/X7T 2 (T7T XL 

T~2 ’ ^ T 2 1 ^ * 

pL^to pLuj z 7r 


From this point onwards, the tildes will be dropped. 

We next proceed to linearize the governing equations for the purpose of making the analysis 
tractable. The typical vertical displacement of the BM is approximately 10 ~ 6 cm which is six 
orders of magnitude smaller than its length of 3.5 cm m This implies that the flow Reynolds 
number is very low, on the order of 10 ~ 6 or less, and so nonlinear effects can be ignored. Assuming 
that there is negligible coupling in the membrane along the longitudinal (x) direction, we only 
consider membrane displacements in the y-direction. Another major simplification is achieved by 
eliminating the delta functions and reformulating the problem in terms of jump conditions across 
the BM. To do so, we integrate the governing equations across the membrane at its linearized rest 
state y = 0 [2SJ, yielding the alternate system of equations 

(2.14) ^ = -Vp+i/Au, 

(2.15) V ■ u = 0, 


away from the BM and 

(2.16) [pj = — Ke~ ax (l + 2 rsinf) h(x,t), 

(2.17) u(x, 0, t) = 0, 

( 2 . 18 ) v(x,0,t) = ^, 

across the membrane where h(x,t) represents the vertical membrane displacement, [p](a;,f) = 
p(x, 0 +,f) — p(x, 0 ”, t) is the jump in pressure across the membrane, and u{x , y, t) and v(x , y, t) are 
the horizontal and vertical components of the vector velocity u. Eqs. (12.141) (12.181) were studied 
analytically in [31] without the time-varying stiffness parameter (that is, with r = 0). 


3. Floquet Analysis. Owing to the presence of a time-varying parameter in the system, we 
invoke Floquet theory [S] to analyze the stability of perturbations of the membrane from its resting 
state. Floquet theory assumes a solution of the form 

(3.1) «(x,f) =e lf P(x,f), 

where the function P(x, t) is periodic in time with period 27r. The Floquet exponent 7 G C 
determines the stability of the solution as t —> 00 . For the purposes of our analysis, we extend the 
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x-domain to [—7r, 7r] and impose even symmetry across x = 0. We then take the unknown variables 
to be of the form 

OO OO 

(3.2) u(x, y, t) = e 7t E E u n k (y) e int e ik *, 

n—— oo k=— oo 
OO oo 

(3-3) v{x,y t t) = e* E E v%(y) e int e ikx , 

n=— oo k=— oo 
oo oo 

(3-4) p(x,y,t) = e^ ]T E pj(l/)e int e ifca! , 

n=— oo k=—oo 
oo oo 

(3.5) = E E h n k e mt e ikx , 

n=— oo k=—oo 

where we have assumed that -P(x, t) can be expanded using a Fourier series in both space and time. 
Our solution approach is similar to that used by Cortez et al. [TO] (modulo the corrections in m) 
who analyzed the stability of a 2D circular membrane in response to a perturbation in the form of a 
single Fourier mode. In contrast, here we must represent each solution mode as an infinite Fourier 
series because of the mode-coupling through the spatial non-uniformity in the stiffness parameter 
which will be shown later in this section. 

We begin by finding solutions for the y-dependent Fourier coefficients and p £. Take the 

divergence of the momentum equations (12.1411 and apply the incompressibility condition (12.151) to 
arrive at a Poisson problem for pressure 

OO 

(3-6) A P = J2 {-k 2 P n k (y)+p n k "(y))£J! = 0. 

n,k——oo 

Here we have simplified notation by setting £ k (x, t) = exp[(y + i n)t + i/ex] and denoting y-deriv- 
atives using primes. The £% are all linearly independent and so we have a decoupled system of 
ordinary differential equations for the p k (y) 

(3.7) ~k 2 Pk(y)+Pk"(y) = 0 Vrq/ce Z, 


which after imposing boundedness in y yields the pressure solution 


(3.8) 


Pk(y) 


n n P k V 

a k e 

h n p~ k V 

°k e 


if y < o, 
if y > o, 


for each n, k G Z, with constants a£ and b £ yet to be determined. Note that taking k = 0 is valid 
in the above expression since this implies simply a constant pressure in each sub-domain. 

We can now solve for the Fourier coefficients of the vertical velocity component v by substituting 
the series (El into the momentum equation to obtain 


(3.9) 


OO OO 

E (7 + i nK(y)£%= E (-Pfc'fo) “ ^v n k (y) + uvfiy))^. 

n,k ——oo n,k =—oo 


This is equivalent to the infinite linear system of ordinary differential equations 
(3-10) vf(y) - m 2 vl(y) = pp n k '(y ) Vn, k G Z, 

where 


Pk 


k 2 


with Re{/3£} > 0. 


(3.H) 
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Assuming that 7 + in 7^ 0 and k ^ 0, the solution is 


2fc/3£ 


(3.12) v%(y) = 


{fik ) 2 - k 2 


n ky , 
a k e + 


k -< + ^rrrk)^ v , if y < 0, 


- k 


i% + k k 


2 vFi 


9 R n 

Zt k h n P ~ k v _ 

m 2 ~k 2 k 


- k 




a k ) e ^ kV i if y > 0. 


The continuity equation reduces to 


(3.13) 


i ku%(y) + vl'{y) = 0 Vra, k € Z, 


from which 11 % (y) are found to be 


2 k 2 


(3.14) ul{y) = 


m 2 ~ k 2 


al e ky 


Pk~k 


a k + 


-6£ I rP* v 


if y < 0, 


2 kv 


2k 


m 2 - k 2 


62 e~ ky + 


k -K + J^rrOe-^, if y > o. 


Pt-k K M + k 


We then impose the continuity conditions on velocity across the membrane 


(3.15) 

(3.16) 


■"(0+) = «2(0") = -2- 


1 


1 


2v V Pk + k 


62 = 0, 


*(0+) = i£(0-) = 


1 


-k 


2vp% \(3% + k K ffi + k 


K = (7 + i«)*fc, 


which can be solved in terms of the as 


(3.17) 

(3.18) 


■m, 


a n k = -f(7 + m) 


ft n -L U 

K = 1/(7 +in) 


k' 


We note that the special cases 7 + in = 0 and k = 0 both yield the trivial solution, ajj = 6£ = 0. 

Using the jump condition (T2.16H . we can now formulate an infinite system of linear equations that 
connects all of the membrane coefficients, h%. Substituting the solution for pressure into (12.161) . 
we obtain 

(3.19) 


E 2^(7 + in) 

n,k=— 00 


on I 7 < - x - ) 

Pj ^—TkK 4”= E -«e-“(l + 2rsin t)h n k £ t 

n,k= — 00 


k 5 


provided that 7 + in ^ 0 and k =£ 0. In the special cases when 7 + in = 0 or k = 0, these equations 
reduce to 


(3.20) 


0= E — ax (1 + 2 t smt)h k £ k . 


,k——c 


In order to proceed any further, we first need to expand the exponential ( e ~ ax ) and sinusoidal (sin t) 
terms in their respective Fourier series. For the time-dependent factor, we can write 1 + 2r sint = 
1 — ire 1 * + ire -1 *, but the Fourier series for the exponential function does not converge uniformly 
because e~ ax is not a periodic function on [0,7r]. It is for this reason that we have extended the 
spatial domain to [—7r, 7r], and then we simply need to remember that is only the portion with 
x > 0 that is of physical interest. We instead use the even periodic extension e~ a ^ of the elastic 
stiffness function on this extended interval. Then, for 7 + in ^ 0 and k ^ 0 we have 


(3.21) 


OO /Q77 I j 

E 2p(7 + in)^-^«£ fc ” 

n,k =—00 


00 

= E -« 

n,k——oo 



(1 — ire 1 * 


+ ire 14 


)K£ 


n 
k > 




























W. KO AND J. M. STOCKIE 


while for 7 + in = 0 or k = 0 , 


(3.22) 


OO 

°= ~ k 
n,k=—oo 



(i- 


ire 14 + ire lt )h 


cn 
k c k > 


where 


(3.23) 


1 - (—l) J e _a7r 
7r (a 2 + j 2 ) 


are the Fourier coefficients of e~ a ' x ' on [— 7 r, 7 r]. The complex exponential terms involving x and t 
in Eqs. (13.211) and (13.221) introduce a shift in the indices of £% in both n and k. which has the effect 
of coupling the corresponding modes. We can then rearrange the sum in order to gather together 
all terms involving the common expression £%, and hence obtain 


9,y2 on 00 00 

(3.24) — (/?fc - k) {ft + kf Bh. K + Y, Cfc-A" = ir Y c ^ h V ~ h T^ 

j= — OO j — ~ OO 


for 7 + in ^ 0 and A: ^ 0 and 

OO OO 

(3.25) £ = ir Y Ck-jih]- 1 - h] +1 ), 


J=-oo 


J=-oo 


for either of the special cases 7 +in = 0 or k = 0. These last two equations comprise an infinite linear 
system for the h™ in which the spatially-dependent stiffness introduces a simultaneous coupling 
between all spatial modes that is not present in the spatially-uniform (a = 0 ) case from | 10 J . 

Because we are interested in investigating the stability of the parametrically forced problem, and 
in particular finding the stability boundary in parameter space, we only have to look for periodic 
solutions of Eqs. (13.241) (13.251) in the situation when Rejy} = 0. Furthermore, there are only two 
distinct values of 7 that are of actual interest for determining the stability boundary: the first 
corresponds to 7 = 0 and will be referred to as the harmonic case; and the second, 7 = ^i, is called 
the subharmonic case. To ensure that the solution h(x, t) is real-valued, we impose reality conditions 
for the Fourier coefficients that apply to both time and space indices. In general, the reality 
condition for a two-dimensional Fourier series is h J? = ft”£ where the overbar denotes the complex 
conjugate. However, we have to consider the harmonic and subharmonic cases separately when 
applying the condition in the temporal modes. Furthermore, we want to ensure an even spatial 
symmetry in our solutions which leads to a decoupling in the reality condition. Consequently, the 


reality conditions that we impose are 


(3.26) 

h- k = K, 


(3.27) 

( h n 

u—n J 11 k > 

k "Ur 1 . 

for the harmonic case (7 = C 
for the subharmonic case (7 


As a result, the reality conditions introduce a certain symmetry between the Fourier coefficients 
such that we only need to consider non-negative values of n and k. 

These choices of harmonic and subharmonic values for 7 can be justified as follows. The solution 
form (13.11) implies that 

(3.28) u(x, t + 2 nn) = e 7 (t+ 2 ™ ) P(x, t) = fw(x, t), 

for any positive integer n where £ = e 727r and t is fixed. As n —> 00 the long term behaviour of the 
solution depends the value of £. We conclude that solutions are stable if |£| < 1 and unstable if 
|£| > 1 , where the special values £ = ±1 correspond to periodic solutions that define the marginal 
stability boundaries separating stable and unstable solutions. If 7 = 0, then £ = 1 and 


(3.29) 


u(x, t + 2n) = u(x, t) 
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which is the 27r-periodic harmonic solution. If 7 = ^i, then £ = — 1 and 

(3.30) u(x, t + 2 tt) = — -u(x, t), w(x, t + 47t) = u(x, t), 


which is the period-doubling subharmonic solution. These same relationships hold for all of the 
dependent variables. 

We are aware of no analytical method for determining solutions to the infinite system (13.241) (13.271) 
explicitly, and so we truncate the system at values ofn = 0,l,...,TV and k = 1, 2,..., M and then 
approximate the solutions numerically. For k = 0, we have that hg = 0 for each n and so there is 
no need to include them in the linear system. The truncated system of equations can therefore be 
represented as a matrix equation 

(3.31) Ah = rBh , 


where 


(3.32) h =[..., Re(K), Im (hi), R e(h n k+1 ), Im(/i£ +1 ), .. ,] T 

is a vector of length 2 M(N + 1) containing all unknown coefficients, A (and B) are block diagonal 
(tridiagonal) matrices respectively, both of block dimension (TV + 1) x (TV + 1) where each block is 
size 2 M x 2 M. The block diagonal matrix A can be expressed as A = diag(7l°, A 1 ,..., A N ) where 
each block has the form 




rc'i, 

+ D? 

Ci ,2 




C 2.1 

C 2 ,2 + D% 


(3.33) 

A n = 







Cm, 1 

Cm ,2 

.. c 

with 






(3.34) 

Ck,j = 


+ Ck-\-j 

0 

0 

Ck—j + e/c+j 

5 

(3.35) 

D n- ^ 
k nk 

' Re{(^ 

. im m 

- km+k) 2 m 
-km+k) a ^} 

For the harmonic case, is simply the 2x2 

zero 


C\,M 

C2.M 




D 


M. 


the form 


(3.36) 


B = 


B B 

B 0 —B 

B 


0 —B 
B 0 


where 


Ci,i 

C\,2 

C\,M 

(3.37) 

B = 

C 2) 1 

C 2 ,2 

C 2; M 



Cm,\ 

Cm,2 ■ 

■ Cm,m. 


Cfc,7 — 


0 

C/c—y T C/c+j 


Ck—j Ck-\-j 

0 


The reality conditions determine the first block row of B. For harmonic solutions, B is the 2 M x 2 M 
zero matrix while B has the same form as B in (13.371) except that the 2x2 sub-blocks are 


Ck,j 


0 2ck~j T 
0 0 


(3.38) 
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For subharmonic solutions B = —B, and B consists of the 2x2 sub-blocks 


Now that the entries in the matrices A and B are all known, we can view (13.31 j> as an eigenvalue 
problem 

(3.40) A~ 1 Bh = —h, 

T 

where the key determinant of stability is the eigenvalue -. In particular, we are only concerned 
with those physically-relevant values of r that are real and are less than 4, since this ensures that 
the elastic stiffness function K(s,t ) is real and non-negative. 

4. Natural Modes for an Unforced Membrane. Before solving the parametrically forced 
problem, we first examine the stability of the unforced membrane corresponding to r = 0. Previous 
stability analyses of the immersed boundary method showed that unforced membranes are always 
stable mm and we expect a similar result here. The Fourier coefficients of the unforced solution 
satisfy 

( 4 ' 1 ) ( k + J~ + fc2 ) h ° + H Ck-jh° = 0 

V \ v / j=-oo 

when k ^ 0 , and 

oo 

(4.2) Y, C *-A° = 0 

j=—oo 

when k = 0. We have introduced the new dimensionless parameter <j> = v 2 /k = n 3 p 2 / (paL 3 ), 
which is a measure of the relative importance of viscous fluid force relative to elastic membrane 
force. The system of equations ED and (14.211 can be written simply as a matrix system Th° = 0, 
where the entries of the matrix T depend on the parameters </>, 7 / 1 / and a. This linear system 
has non-trivial solutions only if det(T) = 0, and so by fixing values of <f> and a we can determine 
values of 7 /^ such that solutions to the homogeneous system exist. In practice, we proceed by first 
truncating the infinite series in nil) and in m to M terms, and then computing the determinant 
numerically. 

Figure |4T] depicts the zero contours of the real and imaginary parts of det(T) for the specific 
choice of parameters a = 1.56 and M = 20, and two values of <p = 5 x 10~ 10 , 2 x 10 -4 . The points 
where the contours intersect correspond to the natural modes of the system. Observe that Rejy} 
is always negative at the intersection points, from which we conclude that all solution modes are 
stable. We see next how the behaviour of the natural modes depends on <p, the relative strength 
of fluid viscosity to membrane stiffness. For the relatively small value of </> = 5 x 10~ 10 , the 
dominant modes (which are slowest to decay) have non-zero Imjy} and are therefore oscillatory. 
When (f> is increased to 2 x 10~ 4 , viscosity has a much stronger influence and the dominant (lowest 
wavenumber) modes decay without oscillations, although decaying oscillatory solutions still do 
exist. In both cases, the unforced modes always decay in time and hence any periodic or unstable 
solutions must arise from a periodic modulation of the membrane stiffness. 

5. Parametrically Forced Pure-Tone Response. In this section, we demonstrate that 
an internally-forced membrane can generate travelling wave solutions that are similar to solutions 
obtained from other models that impose an external forcing. We compute harmonic solutions to 
(13.2411 (13.2511 numerically and then choose the smallest physically-allowable value of r along with 
its corresponding eigenvector. We can reconstruct a periodic solution for h(x,t ) using the series 
representation m , which is then compared to the travelling wave solutions from m obtained 
when a pure-tone external forcing is applied. We take parameter values from )51j (listed in Ta¬ 
ble 15.11) that are also consistent with parameters reported elsewhere in the literature for the human 
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(a) <j) = 5e-10 (b) (|) = 0.0002 




Re{y)/K 1/2 Re{y}/K 1/2 

Fig. 4.1. Zero contours of the real (thick, blue) and imaginary parts (thin, red) of det(T) for: (a) (f> = 5xl0 —10 
and (b) (j) = 2 x 10“ 4 . In both cases, we take a = 1.56 and use M = 20 modes. The resonant modes correspond to 
intersection points between the real and imaginary contours. 

Table 5.1 

Parameter values (or ranges) used in the cochlea simulations. Values are taken from ms and (with the 
exception of a) correspond to the human cochlea. 


Physical parameters (cgs units) 


fluid density 

p = 1.0 gem -3 

fluid viscosity 

/i = 0.02 g cm -1 s —1 

elastic stiffness 

a = 6 x 10 s gem -2 s -2 

elastic stiffness decay rate 

A = 1.4cm -1 

basilar membrane length 

L = 3.5 cm 

forcing frequency 

UJ G [400, 1600] s -1 

velocity scale 

U c = e [ 446 ,1783] ems -1 

pressure scale 

P c - p “ 2 f 2 G [2.0 x 10 s , 3.2 x 10 6 ] dyn cm -2 

Dimensionless parameters: 


dimensionless decay rate 

a = 1.56 

forcing amplitude 

r G [0,0.5] 

dimensionless viscosity 

v G [8.06 x 10 -6 , 1.61 X 10 -4 ] 

dimensionless stiffness 

k G [0.135, 53.9] 

viscosity/stiffness ratio 

d>= = iV j- 4 8 x 10 -10 

^ k pcrL 


cochlea, with the only exception being the elastic stiffness value. Other cochlea models take values 
of a that range from 1 x 10' |2Sj to 2 x 10 9 gcm -2 s -2 pQ; however, we have chosen the smaller 
value of a = 6 x 10 s gcm -2 s -2 for the purposes of this paper in order to permit a direct comparison 
with the results in m for a similar IB model. 

Figure EH] shows solution profiles of the BM displacement curve h(x,t pea k ) for forcing frequencies 
uj = 400, 800, 1200 and 1600 s -1 , where t pea k represents the time when the maximum vertical 
BM displacement occurs. The wave envelope is determined by computing the absolute value 
of a complex-valued function whose real part is the BM profile and the imaginary part is its 
Hilbert transform EHMj. The envelope is normalized so that the maximum height is 1. These 
wave envelopes have the characteristic asymmetric shape seen in experiments |49j . exhibiting an 
amplitude that increases gradually from base to peak, followed by a sharp decline at the apical end. 
The solid curve in each case corresponds to the harmonic mode for which the response frequency is 
equal to the forcing frequency. A qualitatively similar result is obtained for subharmonic solutions 
(dashed curves) except that the response occurs at a frequency equal to half that of the internal 
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forcing and so the wave profiles are shifted. This indicates that the location of the envelope peak 
depends on the response frequency and not the forcing frequency. 

Through a combination of experiments and analysis, von Bekesy showed [54] that the peak loca¬ 
tion of the travelling wave envelope actually depends logarithmically on the stimulus frequency. 
Figure [5^] depicts the computed peak locations at various response frequencies, and our results lie 
nearly on a straight line which is consistent with the logarithmic dependence just mentioned. This 
plot also includes the asymptotic results derived from the model in [31 as a solid line, which is 
clearly very close to our own results. As a further validation, Figure HT3l compares our BM displace¬ 
ment curve with the corresponding result from [31 j for the frequency w = 400 s . Both profiles are 
shifted so that the wave envelope peak occurs at x = 0, after which we can see that the qualitative 
shape and solution envelope are quite similar. These numerical results demonstrate that external 
forcing is not required to obtain a realistic travelling wave solution on the BM, and indeed that 
parametric (internal) forcing in the BM stiffness can generate pure-tone solutions consistent with 
that observed in other models. 


(a) BM profile for to = 400 s’ 



x (cm) 

(c) BM profile for a = 1200 s' 



(b) BM profile for to = 800 s’ 




x (cm) 


Fig. 5.1. Normalized BM displacement profiles for the harmonic (solid) and subharmonic (dashed) cases at 
frequencies ui = 400 , 800, 1200 and 1600 s ~ 1 . 



Fig. 5.2. BM envelope peak location plotted against response frequency, with the horizontal axis (frequency) 
plotted on a log scale. Harmonic (x) and subharmonic solutions (o) have response frequencies ui and ui/2, respec¬ 
tively. 


6. Parametric Resonance. 
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(a) Current Model 



Relative distance from peak (cm) 


Fig. 5.3. Normalized BM displacement profiles from 
(right) at a frequency of lj = 400 s —1 . 


(b) LeVeque et al. 1988 



Relative distance from peak (cm) 

our IB model (left) and model by LeVeque et al. m 


6.1. Case a = 0. To gain more insight into solutions of the eigenvalue problem (13.401) . we 
begin by considering the simple case a = 0 where the elastic stiffness does not depend on BM 
location and thus the spatial Fourier modes are decoupled. For each spatial wavenumber k we have 

(6.1) - k)(J% + + K = i r(h ”- 1 - h n k +1 ) 

for n = 0,l,...,lV. This equation may be rewritten in matrix form as 

(6.2) A k h k = TB k h k , 


where for each value of k 

(6.3) h k = [R e(h° k ), Im(/i£), Re(hf), Im(7if)] T , 

(6.4) Ak = diag(7 + D k , I + ...,/ + D k ), 

D k is defined by (13.351) and I is the 2x2 identity matrix. For the matrix B 1 we need to consider 
separately the two cases corresponding to harmonic solutions where 



' 0 

0 

0 

2 






0 

0 

0 

0 






0 

-1 

0 

0 

0 

1 



B k = 

1 

0 

0 

0 

-1 

0 





0 

-1 

0 

0 

0 1 





1 

0 

0 

0 

-1 0 









and subharmonic solutions where 


( 6 . 6 ) 


0 1 

1 0 

0 1 

-1 0 




0 -1 

1 0 

0 0 

0 0 

0 1 

-1 0 




0 -1 

1 0 

0 0 

0 0 

0 1 

-1 0 








Figure 16.11 contains Ince-Strutt diagrams jB] that depict r (corresponding to the critical forcing 
amplitude) plotted against the spatial wavenumber k for two values of the frequency oj = 500 and 
1000 s -1 . Each point on the diagram represents a periodic solution to the linearized IB equations 
so that when taken together these points trace out the marginal stability contours that divide 
parameter space into regions corresponding to stable or unstable solutions. The contours have a 
characteristic “tongue-” or “finger-like” shape and alternate between harmonic and subharmonic 
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(a) co = 500, K = 2.15, V = 3.22e-05 



spatial mode, k 


(b) co = 1000, k = 0.539, v= 1.61e-05 



spatial mode, k 


Fig. 6.1. Ince-Strutt diagrams showing the critical forcing amplitude r for a = 0 plotted against spatial 
wavenumber when u = 500 s — 1 (left) and u) = 1000 s -1 (right). The physically relevant parameters for instability 
are highlighted by the vertical solid lines. 


solutions, labelled “H” and “S” respectively. Unstable solutions correspond to parameter values 
that lie above and inside each tongue, while parameter values below the tongues yield stable solu¬ 
tions. Although we can expect parametric resonance to occur for any choice of parameters located 
inside one of the unstable tongues, the parameters are further constrained by the following two 
physical considerations: first, that only modes corresponding to integer values of k are physically 
realizable; and second, that the forcing amplitude must satisfy r < i so that the BM stiffness 
K (s, t) remains non-negative. To get a clearer idea of the unstable modes that correspond to phys¬ 
ical BM oscillations, we fix k at integer values ranging from 1 through 6 and display in Figure ROl 
the stability plots as a function of frequency a; on the horizontal axis. In each case the right-most 
tongue is labelled “H” or “S” and the subsequent tongues moving to the left alternate between 
harmonic and subharmonic modes. We observe that as k increases, the contours tend to widen 
and shift downward and to the right; consequently, it is higher frequency modes that are more 
susceptible to resonant instabilities. Indeed, the lowest wavenumber modes can be excited at very 
low values of the forcing amplitude r provided the forcing frequency is high enough. 

To verify the existence of these resonant solutions from our linear analysis in the case a = 0, 
we next perform numerical simulations of the full governing IB equations (12.11) (12.51) . that now 
include nonlinearities from both the advection terms and the delta function integral terms. We 
use a standard approach similar to the algorithm described in [IB] , in which the fluid variables are 
discretized on an equally-spaced rectangular grid and the membrane on a moving set of Lagrangian 
points. A split-step projection method is used to solve the incompressible Navier-Stokes equations 
and a regularized delta function is used to approximate the integral terms that encapsulate the 
fluid-structure interaction. We use the IB algorithm implemented in the freely-available Matlab 
software package Mat IB whose implementation details can be found in [T5] . 

Simulations are performed on a doubly-periodic fluid domain of size [—L, L\ x [— L, L] and we 
use the forcing parameters oj = 900,1000,1100 s -1 and r = 0.1,0.2, with all other parameters 
listed in Table l5Jl Figure [tUTl depicts the time evolution of the peak BM amplitude for an initial 
membrane displacement corresponding to a k = 1 cosine wave with amplitude 10 -6 cm. For the 
parameter values chosen, the results exhibit a range of behaviour including stable (non-resonant) 
solutions in which the amplitude decays over time, as well as resonant solutions that experience 
growth in amplitude by up to two orders of magnitude for the largest values of forcing frequency 
u> and amplitude r. The Ince-Strutt diagram corresponding to k = 1 from Figure 16.21 may be 
used to predict the solution stability in these simulations, and the expected solution behaviour is 
summarized in Table ITH1 for the parameter values corresponding to the simulations. Clearly, the 
linear analysis matches the stability behaviour observed in simulations. Most notably, for the case 
of amplitude r = 0.1 we capture for increasing w how the k = 1 mode transitions from stable, 
through the marginal stability boundary into an unstable tongue, and then returns again to the 
stable region. 
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(a) k = 1 (b) k = 2 








Fig. 6.2. Ince-Strutt diagrams showing the critical forcing amplitude r for a = 0 plotted against forcing 
frequency when k varies from 1 to 6. 


CO = 900, X = 0.1 




CO = 1000, x = 0.1 (resonant) 



\ A A A / 


vvvy 


0 0.05 0.1 


CO = 1000, x = 0.2 (resonant) 



co= 1100, x = 0.1 




Fig. 6.3. Time evolution of the IB peak amplitude when a = 0 given an initial k = 1 mode cosine wave profile, 
for various values of parameters lo and r. 
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Table 6.1 

Analytical stability behaviour predicted for the parameters used in Figure 1 6. 3\ with a = 0. 



ui = 900s -1 

UJ = 1000 s- 1 

UJ = 1100s -1 

T = 0.1 

stable 

unstable 

stable 

r = 0.2 

stable 

unstable 

unstable 


6.2. Case a ^ 0. We next investigate the stability of solutions in the spatially-coupled case 
(a ^ 0) where the BM stiffness varies exponentially along its length. The stability contours are 
shown as plots of r versus w in Figure 16.41 using the same parameters listed in Table 15.11 Here 
we have displayed the harmonic and subharmonic mode plots separately, and we also present two 
sets of results that truncate the series solutions at different numbers of spatial modes, M = 5, 
10 and 250 (in all cases using N = 20 temporal modes). In contrast with the a = 0 results 
from the previous section where stability contours are disjoint and the behaviour of a given mode 
is easy to identify, we observe that contours overlap due to the coupling between spatial modes. 
A similar “mode-mixing” effect has been observed in other physical systems such as the double 
pendulum (2 23. Furthermore, we find that the number of stability contours depends strongly on the 
number of spatial modes M included in the truncated series expansion; in particular, increasing 
M gives rise to more stability contours that tend to pack more closely together. As M gets large, 
the contour “tongue-tips” sweep out a smooth curve that divides parameter space into stable and 
unstable solutions as seen in the bottom row of plots in Figure 16.41 for M = 250. We note that 
convergence in the time modes is much faster than in the spatial modes, thereby requiring that 
M be taken significantly larger than N in order to achieve accurate results. We note in particular 
that when either k > 250 or n > 20, all coefficients satisfy \h^\ < 10 -4 and so the neglected modes 
have negligible effect on our computed results. 

Although it is no longer possible to predict the growth or decay of a single wavenumber mode as 
in the a = 0 case, we can nevertheless still identify the region in parameter space corresponding 
to stable solutions. Figures 16.4b and 16.4F show that for each forcing frequency there is a value of 
r 0 in the interval [0, for which a forcing amplitude 0 < r < t 0 yields a stable solution whereas 
amplitudes r D < r < i lead to resonance. Again, viscosity acts as a stabilizing mechanism in 
the sense that increasing ft will increase r D and consequently increase the size of the region in 
parameter space where solutions are stable. Varying a has the effect of changing the range of 
resonant frequencies: increasing a causes the contours to spread out, thereby increase the range of 
frequencies that result in unstable solutions. 

Once again, we use numerical simulations of the full IB model equations on the domain Q = 
[—L,L\ x [— L,L\ to validate the existence of resonant modes found analytically. For the initial 
membrane configuration, we use the same cosine wave with amplitude 10~ 6 cm and wavenumber 
k = 1. Figures 16.51 16.61 and 16.71 display the time variation of the amplitude of the first three 
Fourier cosine modes with forcing frequencies u> = 400, 600 and 800 s _1 respectively. In each case 
we also choose three different values of the forcing amplitude, r = 0.05, 0.08 and 0.1. Even though 
the initial condition contains a pure wavenumber k = 1 mode, all k- modes are eventually excited 
because of the mode-coupling that arises through the spatially-dependent stiffness. According to 
FigureHHl we expect the r = 0.05 cases to be stable (since the tips of all tongues lie above this value 
of t) while taking r any larger should destabilize the solution. Indeed, numerical simulations with 
r = 0.05 do show that BM oscillations decay in time and that the parametric forcing is insufficient 
to initiate an instability. Furthermore, when r is increased to 0.08, the solutions become unstable, 
sustained oscillations appear, and for the largest value of r = 0.1 the peak amplitude grows even 
larger. It is important to note that in all of the resonant cases simulated, the oscillation frequency 
is half of the internal forcing frequency, which is a common signature of parametric resonance. 

Slight differences arise from the fact that our numerical simulations are on a doubly-periodic domain 
of finite length, whereas the analysis assumes a fluid domain of infinite extent in y. Although we 
have chosen the domain size to be large enough that boundary effects are kept to a minimum, 
there are still interactions between periodic BM copies that cannot be completely eliminated in 
our simulations. 






COCHLEA MODEL WITH PARAMETRIC FORCING 


17 


(a) harmonic, M = 5 
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(c) harmonic, M = 10 



03 

(e) harmonic, M = 250 



03 


(b) subharmonic, M = 5 



03 

(d) subharmonic, M = 10 



03 

(f) subharmonic, M = 250 



Fig. 6.4. Ince-Strutt diagrams for the case a ^ 0, depicting the convergence of the critical forcing amplitude 
as the number of spatial modes M is increased from 5 (top) to 10 (middle) to 250 (bottom). The eigenvalues are 
separated into harmonic (left.) and subharmonic (right) modes. Note the different, axes for the M = 250 case. 


We conclude by discussing the existence of parametric resonances in our cochlea model for the 
full range of physically-relevant BM parameters (namely, er). Figure RT51 displays contours of the 
smallest r resulting in resonance predicted by our analysis using N = 20, M = 250 for a human 
cochlea (left) and a gerbil cochlea (right). The parameters for the human cochlea are taken from 
Table 15.11 except for the membrane stiffness. Experimental values for a reported in the literature 
for human cochleas exhibit a large variation, ranging from a = 1 x 10 7 gem -2 s~ 2 from [29] to as 
high as 2 x 10 9 gcm~ 2 s~ 2 in other two-dimensional models (see [401 Table 1], for example). We 
consider this entire range on the vertical axis of Figure RT51 while the horizontal axis extent covers 
to the the entire range of audible frequencies oj £ [50, 20000] Hz for humans. Stiffness values for 
the gerbil cochlea were extracted from [451 Fig. 8] where several experiments are summarized and 
show a £ [1.5 x 10 8 , 2.0 x 10 9 ] gcm~ 2 s~ 2 (a wider range is plotted). The stiffness decay parameter 
for the gerbil BM is estimated to be A = 3.7cm _1 and the length is L = 1.3 cm [38|. The horizontal 
axis of the right plot covers the audible range for a gerbil u> £ [100, 50000] Hz. From these plots, we 
observe that parametric resonance is possible (r < A) for nearly all parameter values corresponding 
to the human cochlea, except at the highest frequencies and the lowest values of a. Therefore, we 
can conclude that parametric resonance arising from fluid-structure interaction effects is possible 
in our cochlea model for most sounds in the human audible range. A similar conclusion can be 
drawn for the gerbil cochlea. A possible experiment to further explore these parameter ranges is 
to study BM response to in vitro electrical stimulation of the cochlea 0133 wherein electrically 
evoked otoacoustic emissions could be modeled as a parametric forcing in the outer hair cells. 
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Fig. 6.5. Time evolution of the amplitude of Fourier cosine coefficients from numerical simulations for internal 
forcing frequency uj = 400 5“ 1 and stiffness forcing amplitude r = 0.05,0.08 and 0.1 (top, middle, bottom). 
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Fig. 6.6. Time evolution of the amplitude of Fourier cosine coefficients from numerical simulations for internal 
forcing frequency lj = 600 s -1 and stiffness forcing amplitude r = 0.05,0.08 and 0.1 (top, middle, bottom). 
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Fig. 6.7. Time evolution of the amplitude of Fourier cosine coefficients from numerical simulations for internal 
forcing frequency u = 800 s —1 and stiffness forcing amplitude r = 0.05,0.08 and 0.1 (top, middle, bottom). 


(a) 


Human cochlea 


(b) 


Gerbil cochlea 




0.5 


frequency (Hz) 


frequency (Hz) 


Fig. 6.8. Contour plots of minimum r required for parametric resonance with physically relevant parameters 
for the human cochlea (left) and the gerbil cochlea (right). The contour corresponding to r 0 = 0.5 is shown as 
a solid curve, and so only those parameters lying in the lower right comer of the diagram (small stiffness, large 
frequency) correspond to stable or non-resonant solutions. 


7. Summary and Conclusions. A immersed boundary model was developed for the basilar 
membrane in the cochlea, for the purpose of investigating the relevance of parametric resonance 
as a novel mechanism for amplification of BM oscillations. Our model captures the fluid-structure 
interaction that occurs between the basilar membrane and the surrounding cochlear fluid. Our 
work is based upon a previous model from [31J . but includes the additional effects of internal 
(parametric) forcing due to variations in the elastic properties of the BM. The prime motivation 
for introducing such a parametric forcing derives from the work of Mammano and Ashmore [32| 
who have uncovered experimental evidence that oscillations of the outer hair cells embedded in the 
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CP can lead to periodic modulation of the BM elastic stiffness. 

We demonstrated first that a parametrically-forced membrane can produce travelling wave solutions 
that are similar to those observed in [31] for a passive BM. A Floquet stability analysis was then 
used to demonstrate the existence of resonant solutions in the linearized IB equations. The results 
were presented as plots of the marginal stability contours in r-w (forcing amplitude-frequency) 
parameter space. For a spatially homogeneous membrane with a constant value of stiffness, the 
stability contours in the corresponding Ince-Strutt diagram are disjoint for each Fourier mode. 
However, for the realistic case of a BM stiffness that varies exponentially along its length, there is 
a mode-mixing effect in which the stability contours overlap in parameter space. We conclude that 
internal forcing through via the BM stiffness at sufficiently large amplitudes can induce parametric 
instability for any frequency in the physiological range of human hearing. These existence of these 
resonances is verified using numerical simulations of a full two-dimensional immersed boundary 
model of the cochlea. 

Our main conclusion is that parametric resonances arising from fluid-structure interactions in the 
cochlea are worthy of further study as a possible contributing factor in the sound amplification 
ability of human and other mammalian hearing systems. One obvious focus for future research 
is to develop a more complete cochlea model that couples the fluid-structure interaction effects 
(giving rise to parametric resonance) along with an existing model for BM mechanical amplifica¬ 
tion [331133 130] gSJ @5] , which would thereby permit a comparison of the relative importance of the 
combined effects. This would also permit us to replace the BM stiffness parameter J33) with a more 
physiologically relevant (non-separable) function in which the spatial dependence is determined by 
an existing mechanical model that has been validated against experiments. Our immersed bound¬ 
ary model is also an ideal framework to investigate effects such as longitudinal coupling within the 
BM 126] 137] or bending-resistant stiffness that some studies claim are important [SUIT]. Finally, 
the active process in the cochlea has been connected by some authors with spontaneous otoacoustic 
emissions [331 for which a study of the potential impact of parametric resonance on amplifying 
such spontaneous oscillations would be a fascinating topic for future investigation. 

Acknowledgments. We would like to express our thanks to Brittany Froese and Jeffrey 
Wiens for the use of their two-dimensional immersed boundary code MatIB [13]. 
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